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Abstract 

We propose a discrete time dynamical system (a map) as phenomenological model of excitable 
and spiking-bursting neurons. The model is a discontinuous two-dimensional map. We find condi- 
tion under which this map has an invariant region on the phase plane, containing chaotic attractor. 
This attractor creates chaotic spiking-bursting oscillations of the model. We also show various 
regimes of other neural activities (subthreshold oscillations, phasic spiking etc.) derived from the 
proposed model. 



The observed types of neural activity are extremely various. A single neuron 
may display different regimes of activity under different neuromodulatory condi- 
tions. A neuron is said to produce excitable mode if a "superthreshold" synaptic 
input evokes a post-synaptic potential in form of single spikes, which is an order 
of magnitude larger than the input amplitude. While a "subthreshold" synaptic 
input evokes post-synaptic potentials of the same order. Under some condi- 
tions a single spike can be generated with arbitrary low frequency, depending 
on the strength of the applied current. It is called spiking regime. An impor- 
tant regime of neural activity is bursting oscillations where clusters of spikes 
occur periodically or chaotically, separated by phases of quiescence. Other im- 
portant observed regimes are phasic spikes and bursts, subtreshold oscillations 
and tonic spiking. Understanding dynamical mechanisms of such activity in bi- 
ological neurons has stimulated the development of models on several levels of 
complexity. To explain biophysical membrane processes in a single cell, it is 
generally used ionic channel-based models The prototype of those models is 
the Hodgkin-Huxley system which was originally introduced in the description 
of the membrane potential dynamics in the giant squid axon. This is a high 
dimensional system of nonlinear partial differential equations. Another class 
of neuron models are the phenomenological models which mimic qualitatively 
some distinctive features of neural activity with few differential equations. For 
example, the leaky integrate-and-fire model, Hindmarsh-Rose and FitzHugh- 
Nagumo model etc. A new important subclass of phenomenological models is 
the map-based systems. Basically such models are designed with the aim of sim- 
ulating collective dynamics of large neuronal networks. The map-based models 
possess at least the same features of Ordinary differential equations (ODE) mod- 
els, and have more simple intrinsic structure offering an advantage in describing 
more complex dynamics. In order to model basic regimes of neural activity 
we design new family of maps that are two-dimensional and based on discrete 
FitzHugh-Nagumo system in which we introduce Heaviside step function. The 
discontinuity line determines the excitation threshold of chaotic spiking-bursting 
oscillations. For some domain of the parameters, we found on phase plane an 
invariant bounded region containing chaotic attractor with spiking-bursting ac- 
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tivity. The interesting fact is that the dynamical mechanism, leading to chaotic 
behavior of our two-dimensional map is induced by one-dimensional Lorenz-like 
map. We demonstrate also that our model can display rich gallery of regimes of 
neural activity such as chaotic spiking, subthreshold oscillations, tonic spiking 
etc. All these modes play important role in the information processing in neural 
systems. 



I. INTRODUCTION 



The nervous system is an extremely complex system [l| comprising nerve cells (or neurons) 
and gial cells. By electrical and chemical synapses of different polarity neurons form a great 
variety large-scale networks. Therefore, modeling of brain's key functional properties is 
associated with study of collective activity of complex neurobiological networks. Dynamical 
modeling approach j^] is effective tool for the analysis of this kind of networks. First of 
all this approach takes building dynamical models of single neurons. On the one hand, 
such models should describe large quantity of various dynamical modes of neural activity 
(excitable, oscillatory, spiking, bursting, etc.). This complexity is associated with the large 
number of voltage-gated ion channels of neurons. It takes employment of complex nonlinear 
dynamical systems given by differential equations. The canonical representative of this 
type of models is Ho dgkin- Huxley system. It describes dynamics of the transport through 
membrane of neuron in detail. On the other hand, to model neural network consisting of 
the large number of interconnected units it is necessary to create simplified models for single 
neuron to avoid problems that are induced by high dimension and nonlinearity. For example, 
one which is commonly used in simulations is integrate-and-fire model 3|. It represents one- 
dimensional nonlinear equation with some threshold rule. That is, if the variable of the 
model crosses a critical value, then it is reset to new value and the neuron is said to have 
fired. To solve the contradiction between the requirements of complexity and simplicity of 
neuron models phenomenological models were introduced. They describe basic properties 
of neuron dynamics, but these models do not take into account the large number of voltage- 
gated ion channels of neurons. As a rule they involve generalized variable which mimic the 
dynamic of some number of ionic currents at the same time. The examples of this type 
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models are FitzHugh-Nagumo Hindmarsh-Rose [5|, Morris-Lecar {g]. They have the 
form of differential equations systems. However, there is another class of phenomenological 
models of the neural activity. These are discrete-time models in form of point maps. In the 
last decade this kind of neural models has attracted much attention For example 

using a map-based approach Rulkov et. al. [11] have studied dynamics of one- and -two 
dimensional large-scale cortical networks. It has been found that such map-based models 
produce spatiotemoral regimes similar to those exhibited by Ho dgkin- Huxley -like models. 

Neuron oscillatory activity can take a variety of forms One of the most interest- 

ing oscillatory regimes is spiking-bursting oscillations regime, which is commonly observed 
in a wide variety of neurons such as hippocampal pyramidal neurons, thalamic neurons, 
pyloric dilator neurons etc. A burst is a series of three or more action potential usually 
riding on a depolarizing wave. It is believed that the bursting oscillations play crucial role 
in informational transmitting and processing in neurons, facilitate secretion of hormones 
and drive a muscle contraction. This oscillation can be regular or chaotic depending on 
the concentration of neuromodulators, currents and other control parameters. Another in- 
teresting oscillatory regime is an oscillation of membrane potential below the excitation 
threshold, so-called subthreshold oscillation. For example, these oscillations with close to 
10 Hz frequency are observed in olivo-cerebellar system providing highl y c oordinated sig- 
nals concerned with the temporal organization of movement execution 13|, [ij] (see more 
discussion in the conclusion). 

n 

The best known spiking-bursting activity model is the Hindmarsh-Rose system [5|. It is 
three-dimensional ODE-based system involving two nonlinear functions. Spiking-bursting 
dynamics of map-based models has recently been investigated by Gazelles et.al [l5|, Rulkov 



nn 



[16|,ll7|, Shilnikov and Rulkov Il8|,ll9(, Tanaka 



A piecewise linear two-dimensional map 



with a fast-slow dynamics was introduced in [l5(]. It was shown that depending on the con- 
nection (diffusively or reciprocally synoptically) , the model demonstrates several modes of 
cooperative dynamics, among them phase synchronization. Two dimensional map is used for 
modeling of spiking-bursting neural behavior of neuron 17, lisl, 13] • This map contains one 
fast and one slow variable. The map is piecewise nonlinear and has two lines of discontinuity 
on the phase plane. Modification of this model is presented in [isl. The further advance- 
ment of Rulkov model is presented in [l9]. A quadratic function has been introduced in the 
model. Using these modifications authors obtained the dynamical regimes of subthreshold 
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oscillation, corresponding to the periodical oscillation of neuron's transmembrane potential 



below the excitation threshold. In 2l| the dynamics of two coupled piece-wise linear one- 
dimensional monostable maps is investigated. The single map is associated with Poincare 
section of the FitzHugh-Nagumo neuron model. It is found that a diffusive coupling leads 
to the appearance of chaotic attractor. The attractor exists in an invariant region of phase 
space bounded by the manifolds of the saddle fixed point and the saddle periodic point. The 
oscillations from the chaotic attractor have a spike-burst shape with anti-phase synchronized 

ving quasi-periodic oscillation for generating the 



spiking. A map-based neuron model invo 
bursting activity has been suggested in [20]. Izhikevich and Hoppenstead have classified [22 1 
map-based one- and two-dimensional models of bursting activity using bifurcation theory. 

Our goal here is to introduce a new map-based model for replication of many basic modes 
of neuron activity. The greater part of our paper deals with regimes that mimic chaotic 
spiking-bursting activity of one real biological neuron. We construct a discontinuous two- 
dimensional map based on well-known one-dimensional Lorenz-type map 23|] and a discrete 



version of the FitzHugh-Nagumo model ^. This is the system of two ODE: 

X = F{x) — y (1) 
y = e{x- J) (2) 

where x is the membrane potential of the neuron and y is the recovery variable describing 
ionic currents, F is a cubic function of x and J is constant stimulus. This model takes 
into account the excitability and regular oscillations of neuron, but not spiking-bursting 
behavior. We shall introduce a discontinuity in the discrete version for this purpose. We 
find conditions under which this two-dimensional map has an invariant region on the phase 
plane, containing chaotic attractor. In addition we show that depending on the values of 
the parameters, our model can produce phasic spiking and subthreshold oscillations also. 

The paper is organized as follows. In Sec. [TTl we describe the map-based model. Then 
in Sec. IIIII we study one- dimensional dynamics in the case when the recovery variable y is 
fixed. In Sec. IIVI we analyze the relaxation two-dimensional dynamics of the model. Then 
in Sec. |V] we find an invariant region bounding the chaotic attractor in the phase plane of 
the model. In Sec. I VI I we observe other modes of neural activity which could be simulated 
by using this model. 
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II. A MODEL FOR BURSTING NEURAL CELL 



Let / : i?^ — > i?^ be a map (x, y) — > (x, y) of the form 

X — X + F(x) — y — (3H{x — d) 

(3) 

y^y + £{x- J) 

where the x-variablc describes the evolution of the membrane potential of the neuron, the 
y - variable describes the dynamics of the outward ionic currents (the so-called recovery 
variable). The functions F{x) and H{x — d) are of the form 

-moX, a X < drain 
F{x) = < mi(x - a), if Jmin <X < Jmax i^) 
^ -mo{x - 1), a X> Jmax 

{1, if X > 
(5) 
0, if x < 

where 

_ ami J _ mo + ami 

mo + mi mo + mi 

The parameter e {e > 0) defines the time scale of recovery variable, the parameter J is a 
constant external stimulus, the parameters /3 (/3 > 0) and d {d > 0) control the threshold 
property of the bursting oscillations. Here we have chosen this linear piece-wise approxima- 
tion of F[x) in order to obtain a simple hyperbolic map for chaotic spiking-bursting activity. 
However, any cubic function can be also used. The map / is discontinuous map and x = d 
is the discontinuity line of /. We consider only those trajectories (orbits) which do not fall 
within a discontinuity set D — (J^o / ^-^^ where L is the union of points of discontinuity of 
/ and its derivative Df. Besides, we assume that mo < 1, then 

det^^ = l + F'{x)+e>0 
d{x,y) 

for any 5 > and the map / is one to one. We restrict consideration of the dynamics of the 
map / to the following parameter region 

< J < d, Jmin < d < Jmax, rUQ < 1. (6) 



Note that under such conditions we have F'{d) > 0. This condition is very important for 
forming chaotic behavior of the map / as we shall see bellow. For convenience, we rewrite 
the map / in the following form 



/ 



/i, if X < d 
/2, if X > d 



where /i, /2 are the maps 



/i : {x,y) {x + F{x) - y, y + e{x - J)) 

/2 : {x, y) ^{x + F{x) - y- (3, y + e{x- J)). 

III. ONE-DIMENSIONAL DYNAMICS OF THE MODEL 



Let us start with the dynamics of the map / when the parameter e = 0. In this case the 
map / is reduced to a one dimensional map: 



X = X + F{x) ~ yo — PH{x — d) = g{x) 



(7) 



where yo is a constant and it plays the role of a new parameter. The map ([7]) can be rewritten 

as 
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» 


= (1 


- mo 


^X — yo, if X < Jmin 
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[x) 


= qx 


-yo 


— ami, if Jmin < X < d 
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[x) 


= qx 


-yo 


- ami - P, if d<x < Jmax 
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[x) 


= (1 


- mo 


)X -yo + mo- P, if X > Jmax, 



(8) 



where g = 1 + mi. 

Let us fix the parameters a, d, mo, mi and consider the dynamics of the map ([H]) in the 
parameter plane (/?, yo)- We restrict our study of the map / to the following parameter 
region 

(3>f3o (9) 

yo > F{Jmax) - P (10) 
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where Pq = F{Jmax) — F{Jmin)- These conditions allow to obtain interesting properties of 
the map ([3]). Let us find the conditions on the parameter values for which the map / acts 
like a Lorenz-type map |23j]. For that we require that (see Fig. [3]) 

\img2{x) < gsiJmax), \imgs{x) > g2{Jrain)- (11) 

X /'d x\d 

It follows from ffTTl) the following condition on the parameter (3: 

(3 < A, (12) 

where 

I3i = min {q{Jmax - d), q{d - Jmin)} ■ 

The inequalities ([9]) and (|T2l) define on the (rf, /3) plane the region (see FigH]). Let us 
take the parameters d and (3 inside the region, and let us consider the {jS, yo) plane. In 
this plane the inequalities ([9]), ffTOl) and f|T2|) are satisfied simultaneously in region Y. In this 
plane the boundary of Y consists of the three lines (Fig. [2]) 

Bo = {{yo,P) ■.p = Po, yo> F{J^,n)} 

B, = {{yo,(3) ■.(3 = (3i, yo > F{Jmax) - A} 

Tl = {(yo, P)-yO = F{Jmax) -P, Po<P<Pl}. 

Consider the dynamics of the map g for {yo, (3) G Y . This region is separated on four 
subregions by the bifurcation lines 

Do = {{yo,P):yo = F{d), Po < P < Pi} 
H = !^{yo,P):yo = F{d)-^P, Po < P < Pi^ 
T2 = {{yo,P) : yo = F(J™n), Po < P < Pi} , 

corresponding to different dynamics of the map g. The line Do coincides with appearance 
of an unstable fixed point x = a + yo/^i through crossing of the discontinuity point x = d. 
Line T2 corresponds to the fold (tangent) bifurcation of the fixed point x = Jmin (see Fig. 
[3]^a,d)). Line H corresponds to the condition 

\imgs{x) =a + yo/mi. 

x\d 

j/n, P) E H there exists a bifurcation corresponding to appearance of ho- 
2J] to the unstable fixed point. The dynamics of the map g corresponding 



Note that for 
mo clinic orbit 
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to subregions Yi{i = 1, ..4) is shown in Fig|3l If {P,yo) E Yi[jY2 the trajectories of the 
map g tend to stable fixed point x = —yo/mo for any initial conditions different from an 
unstable fixed point (Fig 13] (a), (b)). If {P,yo) G V3lJ^4 ^^e map / has invariant interval 
I = {x : b < X < c}, where 

b = qd — yo — ami — P (13) 
c = qd — yo — ami. 

For parameters (/3, ?/o) ^ ^3 the map g exhibits bistable property, that is there exists two 
attractors, one is a stable fixed point and the second is an invariant set of the interval / whose 
basins of attraction are separated by an unstable fixed point (Fig|3I^c)). For (/3,?/o) ^ ^4 
there exists the interval / (Fig 12] (d)) which attract all trajectories of the map g. 

Check that the map g on the / acts like a Lorenz-type map. The map g will be a 



Lorenz-type if 



23| 



(i) the derivative g'{x) > for any a; G / \ {d}; 

(ii) the set of preimages of the point of discontinuity, D = [J g~^{d) is dense in /; 

n>0 

(iii) \im.g{x) = b, lim g{x) = c. 

x\d X yd 



One can see that (ji]) and (jm]) are satisfied. According to [23| the condition (jnl) is satisfied if 



g'{x)>q>l, xel\{d}. (14) 

For the map g on the interval / we have q = 1 + mi and inequality flT^ is obviously satisfied. 
Therefore the map g on the interval / acts like a Lorenz-type map. The possible structure 
of the invariant set of interval / is controlled by value q. 



Let us find conditions under which the map g is strongly transitive. Recall [23| that a 
Lorenz-type map g is strongly transitive if for any subinterval Jq C / \ {d} there is A; > 
such that 

k 

[JfloDlntl. 

Under the condition flT4l) the sufficient condition for strong transitivity on the interval / are 



23|) 

mm 



in{g"^+\g"2+^} > 2 (15) 
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where ni,n2 G Z+ are such that they satisfy the following conditions 

g,{b),...,g^^{b) G [b,d),g^^^\b) G {d,c] (16) 

g,{c),...,gr{c) G {d,c], gr'-'ic) G [b,d). (17) 

Now let us find condition for the parameter values of the map g under which rii = n2 = 1. 
Consider the condition f|T6l) . Let us take rii = k, where k = 1,2, . . .. It is clear that f|T6|) 
holds if the parameter yo satisfies the following conditions 

yo>Fid)-Pq'{q-l)/{q''+'-l) 

(18) 

yo<F{d)-f3q'^+\q-l)/{q'^+^-l). 
Let us require that (fT8|) for A; = 1 is satisfied for 

{P,d)eB+, 2/0 G Fa U (19) 

From inequalities 0, ffT^ and the definitions of the boundaries Ti and iJ, it follows that 
this requirement holds if 

9 > (20) 

Similarly, for n2 = k we get 

y,>Fid)-Piq'^+'-l)/iq^+^-l) ^^^^ 

yo<F{d)-P{q^-l)/{q^+^~l). 

By the same argument as indicated above we obtain that for n2 = 1 inequalities fl2Tl) hold 
if the conditions ( JT71) are satisfied. For example, let us fix g = 1,65, that is mi = 0.65. In 
this case the map g\i is strongly transitive and therefore it follows from the theorem 3.1.1. 
of |23| that the periodic points are dense in /. We note that all of these periodic points are 
unstable (g > 1) and / is a chaotic attractor. Figl3] (c),(d) illustrates the dynamics of the 
map g on the interval / for regions Is and 14 respectively. 

IV. RELAXATION TWO-DIMENSIONAL DYNAMICS OF THE MODEL AND 
SPIKING-BURSTING OSCILLATIONS 



In this section consider the case e << 1 and J > Jmin- This case corresponds to instability 
of the unique fixed point 0{x = J,y = F{J). Since parameter e is sufficiently small, the 
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dynamic of the map / is a relaxation [25|] similarly by to the case of ODE ([Tj) . The distinctive 
characteristic of these systems is two time and velocity scales, so-called "fast" and "slow" 
motions. Basically fast motions are provided by "frozen" system in which slow variables 
are regarded as a parameters, and it is assumed that small parameter of the system equals 
to zero. Slow motions with size of order of the small parameter are given by evolution of 
"frozen" variable. In case of the map /, x is the fast variable and y is the slow one. Let us 
study the fast and slow motions in our system. 

A. Fast and slow motions 

The fast motions of the model ([3]) is approximately described by the map ([7]). As indicated 
above, the dynamics of the map ([7]) can be both, regular and chaotic according to the 
parameter value (Figl3]). Consider now under conditions (Q, f|T2|) slow motions of the map 
/ on the phase plane (x, y) in the region separated by the following inequalities 

X < Jmini y ^ F{.J'max) (2^) 

In the case e « 1 the motions of the map / have slow features within thin layer Mf (e) 



(thickness is of the order e", < a < 1) [25|] near invariant line 

= {{^^ y) ■ y = - bo, X < Jmin] , 

where 

_ mo f^o ~ , _ sJil-mo + ko) 

ko-— + ^—-e, 6o- . (23) 

Directly from the map / it can be obtained that Wf{e) is invariant line not only for e — > 
but for e < uiq/A also. One can see that the dynamics on the line Wf{e) is defined by 
one-dimensional linear map 

X = (1 - mo + ko)x + bo. (24) 

It is clear that the map has stable fixed point x = J. Therefore for J > Jmin the 
trajectories on Wf{e) with initial conditions x < Jmin moves to the line x = Jmin- AH 
trajectories from layer Mf (e) behave in the same way. 

Let us now consider the stability of the slow motions from Mf (e) relatively to the fast 
ones. Since in the case e = each point of the W^{0) is stable fixed point of the fast map 
([7]) then invariant curve Wf{e) is stable with respect to the fast motions. 
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B. Relaxation chaotic dynamics 



It is follows from the previous description that when e is small enough, the structure of the 
partition of the phase (x, y)-plane into trajectories doesn't significantly change with respect 
to case of equations ([7]), fl24l) . The trajectories of the map / are close to the trajectories of 
( I24|) within the layer of the slow motions near W^{e) and to the trajectories of ([7]) outside 
these layer. Therefore, the motions of the map / are also formed by the slow-fast trajectories. 

Let the initial conditions of / belong to neighborhood Mf[e) of the invariant curve 
Wi{£). Any of these trajectories moves within the layer of the slow motions down to the 
neighborhood of the critical point C : x ~ Jmin,y ~ ^^{Jmin), and continue their motions 
according to the fast motions (see Fig. [3](a)), along y ^ yo = F{Jmin)- Since yo E Y4 the 
trajectory of the map g with initial condition x ~ Jmm tends to invariant interval / (See 
Fig. [3] (d)). Therefore the fast motions of the map / with initial conditions C falls into 
some region D^{e) (see Fig HI), D'^{e) -D^(O) if £ 0, where -D(O) is the parallelogram: 

D^{0) = {{x,y) : 
qd — y — ami — [3 < x < qd — y — ami, (25) 
F{Jma.) -P< y < F{d) - P{q - l)/q} 

In other words the region -D^(O) is one parametrical family yo - indexed of invariant intervals. 
As / is attractor, then D^i^e) is also two-dimensional attracting region. Since the map g 
has interval / for ?/o ^ ^3 U then a trajectory involving the map / belongs to the region 
D'^{e) as long as its variable y do not culminate approximately to the value corresponding 
to the line Hq (Fig. [2]). At the same time variable y is slowly increasing for {x,y) G D~^{£) 
as D^{e) G {x > J}. Thus, within the region D'^{e) the variable y continues to increase 
and variable x evolution is close to chaotic trajectory of the map g. 

Over line H{yQ G Y2) the map g has stable fixed point which attracts all trajectories (see 
FigE] (b)). Hence if the magnitude of the variable y becomes about H then trajectory of 
the map / returns into neighborhood of M^{e). Then the process is repeated. As a result 
of these slow-fast motions the attractor A of the system / {x,y) phase plane appears as in 
(FigUa)). 

To characterize the complexity of the attractor A we calculated numerically its fractal 
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dimension df{A). At appears that df{A) takes non-integer values between 1.35 and 1.9 
(FiglH^b)). Therefore A is chaotic attractor. For the parameter values from Fig. H] (b) 
maximum of the fractal dimension df{A) = 1.8287 is accomplished then J = 0.2661. 

V. INVARIANT REGION, CHAOTIC ATTRACTOR AND SPIKING- 

BURSTING OSCILLATIONS 



Let us prove that the system ([3]) has an attractor A for different values e and let us find 
conditions under which the map / has an invariant region. To do that, we construct some 
ring-like region S. Denote by F the outer boundary and by 7 the inner boundary of the S. 
The S is an invariant region if from the conditions (x, y) & S and (x, y) ^ D follows that 
{x, y) G 5*. Its should be fulfilled if 

(i) the vector field of the map / at the boundary F and 7 is oriented inwards to S\ 

(ii) the images /i(F), /i(7), fi{d),i = 1, 2 of the boundaries F, 7 and the discontinuity line 
d belong to S. 

We construct boundaries F and 7 in the form of some polygons. Taking into account the 
condition ([1]) and analyzing the vector field of the map / at the lines with uncertain slope we 
have found the shape of F and 7 (see Fig. [6] (a)). The equations of the boundaries of F and 7 
are presented in the Appendix. Analysis of the position of the images /i(F), /i(7), fi{d), {i = 
1,2) on the phase plane {x,y) show that the condition (In]) holds if {P,d) G (see section 
mil) and inequalities 

(J - Jmin) - - ki{d - J) > F{J^ax) " 

mo 



\/e < min 



B mo[B ~ mi{Jmax - d)] 



2{d-jy J-Jrr 



\fE B B'^ B 
— (J - Jmin) + mi{Jmax - d) + mi{Jmax ^ d) < — + \ — e{d - Jf , if i/i < 



\" "mini I " " L\" max i " " max ^ n ' \ a ")i'-'-v^^r\/7 t\ 

mo 2 V 4 2(a — J) 
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- P - rni{Jmax - d) / T / (^(l +mi) - ami - [3 

^{Jmax -d) I 

(l + mi)(J- J^,„) ^ + (1 + (rf _ J) _ ^ > 0, 

d + F{d) + mo J - Jmax + —{d - Jmin) + ki{d - J) < 

mo 



-y/e < min 



mo mi mi((i — J) 



2, 2, J Jmin 

are satisfied (the parameters fci, 5 and Jo have been introduced in Appendix). Figl6]^b),(c) 
illustrates the transformation of S by the action of the map / under conditions (!26l) . The 
inequahty fl2B]) determine the parameter region Di^v in the parameter plane (J, e) (Fig. [7] 
(a)). Since 

y<y, for (x, y) G 5 P|{x < J] (27) 
y>l/, for (x,i/)G^f|{x> J} (28) 

then the trajectories with initial conditions (x, y) (z S execute rotation motion around the 
fixed point O forming some attractors A. We calculated numerically fractal dimensional 
df{A) (FigJ7|(b)) in terms of e. Its shows that A is chaotic attractor. The possible structure 
of the attractor A in the phase plane is shown in Figl8](a). Fig. [8](b) illustrates time evolution 
of the variable x corresponding to chaotic attractor A. It shows chaotic spiking-bursting 
neural activity. Fig. [7](b) shows that fractal dimension df{A), on average, tends to decrease 
with increasing e. There is a critical value, e = 0.0461, for which fractal dimension has a 
minimal value df{A) = 1.5114. The mechanism of this decreasing can be accounted for by 
the different types of the dynamics of the variable y for different e. As the parameter e 
increases, the velocity of the variable y is expected to climb. Therefore "life time" of the 
trajectories in the strip corresponding to Lorenz-map dynamics is reduced. As a result, the 
chaotical motions are reduced. 
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VI. THE GALLERY OF THE OTHER ATTRACTORS AND THE REGIMES 
OF NEURAL ACTIVITY 



At previous sections it was shown that system ([3]) allows to simulate spiking-bursting 
behavior of the neuron. Here we show that other regimes of the neural activity (phasic 
spiking and burstings threshold excitation, subthreshold oscillation, tonic spiking and chaotic 
spike generation) can be obtained by using the map / also. To do that we neglect the first 
inequahty in ([6]), inequality ([9]) and condition y > F{Jmax) — P- 



A. The generation of phasic spikes and bursting 

Studying response of the neurons to the influence of external stimulus is one of important 
task of neuroscience l| associated with the problem of information transmission in neural 
system. Usually external stimulus is represented as the injection of electrical current into the 
neuron. Let us suppose that the neuron is not excited initially, that is, it is in steady state 
(rest). In the model ([3]) such state of neuron corresponds to stable fixed point O. Consider 
the response of the system ([3]) to pulse type stimulus. We assume that the duration of each 
pulse is small enough (see Fig. [9](a)) and its action is equal to the instantaneous changing of 
the variable x on the pulse amplitude. Besides, we suppose here that e << 1 and therefore 
the dynamics of the system ([3]) is a relaxation. For this parameter region the system ([3]) 
has two thresholds. The first threshold is determined (see Fig. M^c)) by the thin layer of the 
slow motions near the following invariant line 

^r(^) = {(X^y) ■ y = klX - hi, Jmin <X <d} 

where 

, eJ , mi /m? 

Analogously, the second threshold is defined (see Fig. MyC)) as the thin layer of the slow 
motions near invariant line 

W2{e) = {{x,y) : y = kix - b2, d < x < Jmax] , 

where 

02 = mia - hp. 

h 
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Denote by the x'^ and y'^ {y'^ ~ F{J)) the values of the variables x and y after stimulation 
respectively. Let E be the trajectory of the system ([3]) with this initial conditions. In other 
words E is response of the system to pulse input. 

(i) If the amplitude of stimulus is not enough (Fig. [nia),(i)) for overcoming the first 
threshold, then the maximum of the response will be about amplitude of the stimulus. 
Therefore, in this case the generation of the actions potential does not take place. 

(ii) Let us increase the amplitude of the stimulus as it breaks the first threshold but at 
the same time it is not enough for overcoming of the second threshold (Fig. [9](c),(ii)). In this 
case the fast motions of the map / will be close to the fast motions of the nap g on interval 
/ for 1/0 ^ ^3 U ^4- And so, the trajectory of the map / perform some number of irregular 
oscillations around discontinuity line x = d (Fig. [n^c)(ii)). After that, the trajectory E 
within layer near Wl{e) tends to fixed point O (Fig. [9](c)(ii)). Such trajectory E forms the 
region of phasic bursting activity 22| with irregular number of spikes (Fig. [9](b),(ii)). 

(iii) If the amplitude of the stimulus (Fig. [9](a),(iii)) is enough for overcoming the second 
threshold, then the point x = x'^,y = y'^ belongs to the region of attractor of the invariant 
line W2{e), where 

(^) = {{X^y) ■ y = -koX -h, X> Jrnax} 

with 

63 = -mofl - — + p. 

ko 

Therefore trajectory E tends to thin layer of slow motions near invariant line W2{s). It 
moves within thin layer to the neighborhood of the point {x ~ Jmax,y = F{Jmax) ~ 
(Fig. [n](c),(iii)) and its motions continue along fast motions. These motions are close to the 
trajectories of the map g for yo G Yi (see Fig. [3]^a)). Therefore the trajectory E tends to 
the layer near stable invariant line WKe). After that the trajectory E moves within thin 
layer near W2{e) and it tends to the fixed point O (Fig. [9]^c),(iii)). In this case trajectory 
E corresponds to phasic spike [3] (Fig. [9]^b),(iii)). 
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B. Oscillatory modes of the neural activity 



1. Close invariant curve and subthreshold oscillations 

Let us consider dynamics of the map / under following conditions on the parameters. 

e < mo, ruQ > m\/A 

f ml rri^ 
e > max < — ^, — 
(44 

One can see from the Jacobian matrix that in this case the fixed point O has a complex- 
conjugate multipliers. This point is stable for J < J^m and unstable for J > Jmin- Therefore 
the piece-wise map / produces Neimark-S acker like bifurcation (in classical case of Neimark- 
Sacker bifurcation the map is smooth). The fixed point O is surrounded for J > Jmin by 
an isolated stable attracting close curve Cth (Fig. [TOl(a)). The oscillations corresponding to 
the Cth occur under the threshold of excitability of the neuron and therefore it is called in 
neuroscience Q, subthreshold oscillations (Fig. [TOT b)). 

2. "Two- channel" chaotic attractor and chaotic spiking oscillations 

Let us consider again the relaxation {e << 1) dynamics of the map / in the case J > Jmin, 
that is when fixed point O is unstable. Additionally we assume that the parameters of the 
map / has satisfied the following conditions 

F{Jmin) > F{d) - /?, (29) 



FiJma.) > F{d). (30) 

In this case the invariant line TV^(e) separate the fast motions on two flows (Fig. [TlT a)) 
in the neighborhood of the discontinuity line x = d. The first flow forms the trajectory 
performing the chaotic oscillations near the line x = d (Fig. [TlT a)). Their dynamics are 
close to the dynamics of the map g on interval / for yo G |J 14. The second flow consists of 
the trajectories overcoming the second threshold (Fig. [TTlfa)). It moves to neighborhood of 
the stable invariant line W2{e). After that these trajectory tends to the stable invariant line 
Wl{e) and the described process is repeated. These trajectories form chaotically switching 
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flow from one to otlier. As a result is the appearance of a chaotic attractor Ath on phase 
plane (Fig. [TlTa)). The fractal dimension df{Ath) = 1.30335. The attractor Ath determines 
chaotic regime of spiking activity over the background of the subthreshold oscillations (Fig. 

nKb)). 

3. Close invariant curve and tonic spiking 

Let the parameters of the map / satisfy the same conditions as in the case of previous 
subsection IVl B T\ with exception of inequality (12^ . In this case the parameter f3 is small 
enough. Therefore the trajectories with initial conditions from neighborhood of W^i 2(^) do 
not change direction of motion when they intersect the discontinuity line x = d. This leads 
to the appearance on the phase plane of the motions between layer near W{{e) and W2{£)- 
Such dynamics leads to forming close invariant curve Cgp (Fig. [TWa)). So there exists only 
one attractor on the phase plane formed by this invariant closed curve. This determines 
tonic spiking regimes of neural activity (Fig. [T2l(b)). 

VII. CONCLUSION 



A new phenomenological model of neural activity is proposed. The model can reproduce 
basic activity modes such as spiking, chaotic spiking-bursting, subthreshold oscillations etc. 
of real biological neurons. The model is a discontinous two-dimensional map based on the 
discrete version of the FitzHugh - Nagumo system and dynamical properties the Lorenz-like 
map. We have shown that the dynamics of our model display both regular and chaotic 
behavior. We have studied the underlying mechanism of the generation of chaotic spiking- 
bursting oscillations. Sufficient condition for existence chaotic attractors in the phase plane 
are obtained. In spite of idealization, the dynamical modes which are demonstrated in 
our model are in agreement with the neural activity regimes experimentally found in real 
biological systems. For example, subthreshold oscillations (see Fig. ITOlfb)) is a basic regime 



network which plays a key role 



of inferior olive (I.O.) neurons 13 1. Inferior olive neurons belong to the olivo- cerebellar 



14( 1 in organization of vertebrate motor control. It is also 



typical for I.O. neuron a spiking regime over the chaotic subthreshold oscillations (see 
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Fig. [TT]) . The spiking-bursting activity is significant for many types of neurons, in particular 
in hippocampal pyramidal cell |27| and thalamic cells 28|. 

The table summarizes results on gallery of behavior of neural activity showed by our 
model. We hope that our model will be useful to understand the mechanism of neural 
pattern formation in large networks. 
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Parameters 



TABLE 



The regimes of neuronal activity 



J ^ J'min ) £^ ^ ^ 1 



F{J)<F{Jraa.)-P, 

F{J^ax) - P > F{d) (spike) 
F{J) > F{d) - 13 (bursts) 



Phasic spikes and bursts 





AAA 









J > Jmini £ < '^0 

mo > ml/4: 



e > max 



I ' 4 J 



Subthreshold oscillations 



Inequalities ([26 



Chaotic bursting 
oscillations 




F{J,mn) > F{d) - (5 

F^Jma.) > F{d) 



Chaotic spiking 



F(J™„) <F(rf)-/3. 



Tonic spiking 
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APPENDIX. THE EQUATIONS OF THE BOUNDARIES OF THE INVARIANT 
REGION 

The boundary F is given by 

{x,y) : d < X < Jmax, y = {J - Jmin) - m^J -ki{d- J) 

mo 

T2= Ux,y) : J <x <d, y = -p{x -J) (J - Jmin) - rrioJ 

I mo 



with 



r3 = 




• J rain ^ X y 


= {X - Jmin) - moJ } 






mo J 


r4 = 


{{x,y) ■■ 


Jl X J min 1 


y = -ko{x - Jmin) - mo J} 




{{x,y) ■■ 


: — Jl < X < —Jo, 


^/ey ^ {x + Jq) + ^/emoJo} 


r6 = 


{{x,y) ■■ 


Jq ^ ^ Jmin) 


y = mo Jo} 




{{x,y) ■■ 


J min ^ ^ ^ J, y 


= -2^/e{x -d) + F{d) + ^/s{d 




{{x,y) ■■ 


X = Jmax, yi<y < 


y2} 



yi = (J - Jmin) - moJ - ki{d - J) 

mo 

y2 = -2VE{Jmax -d) + F{d) + V^{d - J) 



p 



^ -f r: ^ B 



— £, if y/e < 



Jo = 

Jl = 



2(d-J) Y 4((i-J)2 ^' '■^ — 2((i-J) 

B-mi{Jmax-d) -r ^ B 

d-J ' n yt ^ 2{d-J) 

B = (3- F{d) - mo J, 

V^[{d-J) + 2{d-Jmin)]+F{d) 



mo 

1 + ^/emo)Jo - \fe{koJmin - moJ) 
1 + ko\fe 
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The boundary 7 is given by 

71 = {{X, y) : Jmin <x< J, y= -moJmin} 

72 = {{X, y)-X^ Jmin, -moJmin <y < mi{J - o)} 

73= {{x,y) : Jmin < X < J, y ^ mi{J - a) + ^/£{x - Jmin)} 

74= l{x,y) : J < X < J + ^^^^^^ ^^mm)^ y ^ -ui^(^J _ a) + y/e{J - Jmin)\ 

I ^1 J 

J / \ 7 . \^\J Jmin) ^ ^ \ 

75 = < (x, 7/) : a; = J H , yj,<y<y^) 

76= \ {x,y) : J <x < J + — ^^mm)^ J/ = _ J) _ ^qJ^.^I 
I mi J 



with 



£(>/ Jmin) J 

VS — ^0<^mm 

mi 

y4 = mi( J - a) + V£(J - J„i„) 
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FIG. 4: (a) The chaotic attractor A on the phase plane {x,y); (b) Waveform of relaxation spik 
bursting oscillations generated by the map /. Parameter values: J = 0.13, mo = 0.4, mi = 0.65, a 
0.2, d = 0.3, /5 = 0.25, e = 0.002. 
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FIG. 5: Fractal dimension df{A) of the attractor A versus parameter J. Parameter values: 
mo = 0.4, mi = 0.65, a = 0.2, d = 0.3, /3 = 0.25, e = 0.002 
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FIG. 6: (a) The ring-like invariant region S. (b) The image of the region S f^{x < d} under action 
of the map /i; (c) The image of the region S'Pljx > d} under action of the the map f2- 
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FIG. 7: (a) Parameter region Dinv on the parameter plane (J, e), parameter values: tjiq = 0.5, rrii = 
0.65, a = 0.2, d = 0.34, (3 = 0.31. (b) Fractal dimension df{A) of the attractor A versus parameter 
e, J = 0.15 
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FIG. 8: (a) Chaotic attractor A on the phase plane (x, y). (b) Spike-bursting oscihations generated 
by the map /. Parameters value mo = 0.5, mi = 0.65, a = 0.2, d = 0.34, /3 = 0.31, J = 0.15, e = 
0.004 
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FIG. 9: Response of the model ([3]) to positive pulse e(n): (a) - three different amplitude of the 
stimulus; (b) - the behavior of variable x (membrane potential of neuron); (c) - the phase plane. 
Parameter values: J = 0.119, e = 0.004, d = 0.25, (3 = 0.19, mo = 0.4, mi = 0.8, a = 0.2 
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FIG. 10: (a) Close invariant curves on the phase plane {x,y). (b) Subthreshold oscillations 
generated by map /. Parameter values: d = 0.3, tjiq = 0.4, mi = 0.3, a = 0.2, J = 0.08572, e = 
0.025, /3 = 0.3 
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FIG. 11: (a) "Two-channels" chaotic attractor Af^- (b) Chaotic spiking against the background 
subthreshold oscillations. Parameter values: J = 0.1123, e = 0.004, d = 0.3, mo = 0.4, mi = 
0.3, a = 0.2,(3 = 0.09 
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FIG. 12: (a) Close invariant curve, (b) Tonic spiking. Parameter values: J = 0.1123, e = 0.004(i = 
0.3, mo = 0.4, mi = 0.3, a = 0.2, /? = 0.05 
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